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Abstract —This paper proposes a method to embed the AC 
power flow problem with voltage magnitude constraints in the 
complex plane. Modeling the action of network controllers that 
regulate the magnitude of voltage phasors is a challenging task 
in the complex plane as it has to preserve the framework of holo- 
morphicity for obtention of these complex variables with fixed 
magnitude. Hence this paper presents a significant step in the 
development of the idea of Holomorphic Embedding Load Flow 
Method (HELM) Q], introduced in 2012, that exploits the theory 
of analytic continuation, especially the monodromy theorem @ 
for resolving issues that have plagued conventional numerical 
methods for decades. This paper also illustrates the indispensable 
role of Pade approximants for analytic continuation of complex 
functions, expressed as power series, beyond the boundary of 
convergence of the series. Later the paper demonstrates the 
superiority of the proposed method over the well-established 
Newton-Raphson as well as the recently developed semidefinite 
and moment relaxation of power flow problems. 

Index Terms —AC power flow, voltage control, power flow 
feasibility, voltage stability, holomorphic functions, monodromy, 
analytic continuation, continued fractions, Pade approximants, 
Stahl’s compact set. 

I. Introduction 

Power flow, the most fundamental concept in power system 
engineering, is at the heart of studies ranging from daily 
operation to long-term planning of electricity networks. AC 
power flow problem is a system of nonlinear algebraic equa¬ 
tions that mathematically models the steady-state relations 
between the phasor representation of parameters and unknown 
states in an AC circuit. The parameters typically consist of 
power generated and consumed by source and sink nodes 
and the electrical properties, i.e. the impedance, of lines that 
connect these nodes. The unknown states are primarily voltage 
phasors but could also include continuous or discrete variables 
associated with network controllers, e.g. FACTS devices and 
tap-changing or phase-shifting transformers. The accurate and 
reliable determination of these states is imperative for control 
and thus for efficient and stable operation of the network. 
In certain studies it is equally vital to determine for which 
parameter values the power flow problem becomes infeasible 
as this condition is intimately linked to saddle-node bifurcation 
and the voltage collapse phenomenon 0-0. The distance in 
the parameter space to power flow infeasibility can serve as a 
margin of voltage stability 0. This is certainly one of the most 
theoretical areas in electrical engineering. As conventional 
power systems undergo a fundamental transformation by large- 
scale highly-variable wind and solar generation, distributed 
across the network, this field can experience a resurgence 0. 


Given the inherent limitations of traditional methods, deeper 
understanding of these complicated phenomena requires new 
theoretical approaches rooted in complex analysis and alge¬ 
braic geometry. 

The basis of power flow is Kirchhoff’s current law which 
states that for every node i in Af, the set of all nodes, /,, the 
net current flowing out of that node, is related to its voltage 
V, and those of its adjacent nodes 14 in the following way: 

h= X ** = X X (1 > 

k£j\f{i) k£Af(i) ^ k£j\f[i] 

Af(i) and Af[i\ are the open and closed neighborhoods of 
node i. In- is the current flow through the line connecting 
node i and k and Z,k = Rik + j^ik is the impedance of that 
line which is used to construct the diagonal and off-diagonal 
elements of the admittance matrix as, 
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Since complex power is Si = Pi + jQi = Vil*, the power 
flow problem in its complex form can be expressed as, 


s*= WkYik VieAf-W (3) 

k^M[i] 

Here r is the voltage reference node with |I4| = constant 
and arg(U r ) = 0. It also serves as the slack node meaning 
that S r is a free parameter that accounts for the mismatch of 
complex power and its losses throughout the network. 

The numerical methods, developed historically to solve this 
problem, take the polynomial system of 0 out of its complex 
form by reformulating it either in rectangular form as in 0 or 
in polar form as in 0 where V) = ei+jfi = \Vi\ exp (jdi) and 
Yik = Gik + jBik . These techniques, all based on Newton’s 
method or its variants, iteratively linearize and approximate 
Pi and Qi in 0 or 0, starting from an initial guess. 
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There are two inherent shortcomings in such methods that 
can arise near the feasibility boundary of © or (J5J character¬ 
ized by the saddle-node bifurcation manifold in their parameter 
spaces. Physically, proximity to the feasibility boundary cor¬ 
responds to a network operating close to its loadability limit 
such as periods of peak electricity demand. The first issue is 
the poor convergence of the tangent-based search methods, 
likely due to the increased curvature of the hypersurfaces 
in © or ©. The second issue near the feasibility boundary 
where two or more algebraic branches coalesce at a branch 
point is convergence to solutions that lie on other algebraic 
branches. Although dependent on the dynamical model of the 
physical system, these solutions in power systems typically 
signify unstable Q or low voltage Q operating points. Most 
of these operating points cannot be physically realized and are 
thus false solutions. The region of initial guesses in Newton’s 
method that converge to a particular solution has a fractal 
boundary. The multiple fractal domains of convergence are 
pressed together near the bifurcation manifold which explains 
erratic behavior of such methods in finding the desirable, i.e. 
stable/high voltage, solution even with seemingly reasonable 
initial guesses ©• 

Recently a semidefinite relaxation of rectangular power flow 
in 0 has been reformulated as a special case of optimal 
power flow where the objective function of the semidefinite 
programming (SDP) is minimizing active power loss ©. This 
addresses the convergence failure of iterative methods but has 
its own serious drawbacks. First, the relaxation may not be 
tight and yield a high rank matrix where it is impossible 
to directly recover any solution to the original power flow, 
let alone the desirable one. Second, if the solution of the 
relaxed problem is high rank, nothing can be concluded on 
the feasibility of the power flow in the same vein as non¬ 
convergence of iterative methods cannot rule out the existence 
of solutions. Third, the suggested heuristic, i.e. active power 
loss minimization, does not always find the stable/high voltage 
solution branch. The first two problems can be remedied, at 
least in theory, by obtaining higher-order and thus tighter 
relaxations of ©. The computation cost, however, explodes 
with the order of relaxation and the number of variables. 
Reference na discusses the theoretical underpinning of this 
approach in the context of the generalized moment problem 
and highlights its connection to real algebraic geometry which, 
we see as an obstacle to distinguishing the desirable solution 
branch for algebraic problems of a complex analytic nature. 

Among the above issues, the challenge of finding the 
solution on the desirable branch, more than anything else, 
underlines the significance of embedding the power flow prob¬ 
lem in the complex plane where the extraordinary potentials 
of analytic continuation theory for multi-valued complex func¬ 
tions can be tapped. This is pioneered by the idea of holomor- 
phic embedding load flow (HELM) which builds on the fact 
that under no load/no generation condition (.S', = 0 Vi G A f), 
the network has a trivial non-zero solution for voltage phasors. 
This corresponds to all currents In- being zero and reference 
voltage V r propagated across the network and defines the germ 
of the stable/high voltage branch. Analytic continuation of 
this germ is guaranteed by monodromy theorem to yield the 
desirable solution all the way to the closest bifurcation in the 
parameter space of © where there is a non-trivial monodromy 
and the physically meaningful solution ceases to exist. 


Although the idea of HELM has aroused significant interest 
in the power system community, it yet has to prove its superi¬ 
ority over conventional methods. In this paper we demonstrate 
how the magnitude of complex variables can be held fixed 
while preserving the framework of holomorphicity. This is an 
important step in the embedding of power flow as it models the 
action of network controllers in the complex plane. We also 
show the indispensable role of Pade approximants especially 
for the case of voltage controllers that fix the magnitude but 
not the argument of voltage phasors. Throughout the paper 
we refer to this method as PA to highlight the central role 
of rational approximation of functions of a complex variable 
for recovering the power flow solution. With this abbreviation 
we also want to emphasize the critical direction of research 
for further development of this method. In Section II we 
succinctly review the main ideas of HELM as presented in the 
original paper ID, i.e. for the PQ buses, introduce the concept 
of rational approximation of analytic functions in relation to 
power series and continued fractions and explore the zero- 
pole structure of Pade approximants for a 3-bus example. In 
Section III we introduce the mathematical static model of 
the most prevalent controller in the network, the automatic 
voltage regulator (AVR) of the generator. We demonstrate 
through modification of the previous 3-bus network, this 
time with a generator (PV bus), how the approximation of 
functions of a single complex variable is essential for analytic 
continuation of the voltage phasors with fixed magnitude. We 
interpret the zero-pole structure of Pade approximants and its 
transformation as the solution reaches the feasibility boundary 
and explain the significance of the zero-pole distribution of 
the Pade approximants in terms of voltage stability margin 
at a given operating point. In section IV we demonstrate the 
superiority of the PA over conventional and recently developed 
methods of solving the power flow problem. We introduce a 
7-bus network where Newton-Raphson either fails to converge 
or converges to an unstable/low voltage solution as the active 
power output of a given generator changes. We also show 
that in this network first-order semidefinite relaxation only 
obtains the solution in a small subset of the stable solution 
branch and the second-order (moment) relaxation obtains the 
false solution branches. The numerical results of these methods 
are contrasted with that of the PA method which consistently 
obtains the stable/high voltage solution whenever it exists and 
declares the non-existence of a physically meaningful solution 
beyond the closest saddle-node bifurcation. In this network 
the zero-pole distribution of the Pade approximants depicts 
the analytic structure of the voltage phasors and confirms the 
general pattern of voltage stability margin observed in section 
III. In section V we summarize the key contributions of the 
paper and outline the ongoing work. 

II. Embedding the System of Equations in the 
Complex Plane 

Consider the following parametrization of © in terms of 
sgC with V* replaced with independent variables W t , 
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From a geometric point of view, the 2 n equations of © de¬ 
fine generically an affine algebraic curve in (s, Vi, V 2 , ..., W n ). 


It follows from the Kirchhoff’s current law and the existence 
of the voltage reference node (with V r appearing in ([6]> as a 
parameter) that the polynomials on the right side of (I6at-(l6bl> 
(i- e - E^i] WiV k Y ik and E fce ^[i] ViW k Y& Vi eJV- {r}) 
are algebraically independent. To establish this algebraic in¬ 
dependence in relation to the reference node requires rigorous 
analysis, a task which lies outside the scope of this paper. 
Taking the algebraic independence of these polynomials for 
granted, degenerate cases where the equations of 0 define 
not an algebraic curve but a higher-dimension algebraic variety 
can only arise when the power flow problem is ill-defined as in 
the case of networks with disconnected graphs. This is in line 
with the physical intuition that in the absence of a reference 
voltage, voltages are floating and a given V t can assume any 
value in C. The equations of ([6l> generate an ideal and give a 
starting basis for finding the corresponding reduced Groebner 
basis DU. For any lexicographic order, such as ... > Vi > s, 
this gives a last basis element which is a bivariate polynomial 
fi(Vi,s) for well-defined problems. fi(Vi,s) = 0 can be 
solved for an algebraic (multi-valued) function V-, = V,(sj 
which has holomorphic branches where dfi(V.,s)/dVi ^ 0. 
By permuting the order, we arrive at 2 n algebraic functions 
Vi = Vi(s), Wj = Wj(s), i = 1 giving an algebraic 

parametrization of the curve defined by (0. A branch point of 
this curve will be where any of the components Vi (s) or Wi(.s) 
has a branch point, i.e., fi(V,, s) = 0 and dfi(V,, s)/dVi = 0 
(and similarly for the Wj). The branch point closest to the 
point so at which the Taylor series expansion of any single¬ 
valued branch is developed, determines the radius of conver¬ 
gence of the series. Branch points play a critical role in the 
analytic continuation of these solutions and the PA method. 
This analysis also extends to the case of voltage control where 
we introduce new variables S)(s) and Vi(s). 

Now that Vi(s) is analytic in s, (Vi(s*))* is also analytic 
in s and identical to the conjugate of Vi(s) on the real axis. 
Hence solving (0 is equivalent to analytic continuation of the 
solution of the following system from s = 0 to s = 1, 

= E V* e M - M (7) 
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By defining V(s) = E^=o c ™ s "’ l / V i( s ) = EE=o 
and l/(V)(s*))* = ErEo EE", this system is adequately 
described by the following set of power series relations. 
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The procedure to obtain the coefficients of ® starts by 
setting s = 0 in (l8al >. This gives the linear system of 
EfeeAffi] Y 'i'kCo [k] — 0 which always yields the trivial solution 
(i.e. Cg = V r Vz G Af — {?’}). Next = l/c$ by setting 
s = 0 in < [8b1) . The higher order coefficients are progressively 
obtained by solving the linear system of < [9ab (Vz G A/” — {?’}) 
which itself is obtained by differentiating dial) with respect to 
s and evaluating at s = 0 and the convolution formula of ( l9bl ). 
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The radius of convergence is R = lim r! _ s . 00 |c„|/|c n +i|, if 
the limit exists. This marks the distance from the origin to the 
closest branch point. Notice that when an analytic function 
does not have a closed-form expression as in this case, its rep¬ 
resentation as a power series expansion can be approximated 
by a partial sum of a finite order. Since this approximation 
for V,(s ) does not converge for |s| > R the analytic con¬ 
tinuation of these complex functions toward s = 1 requires 
an alternative representation of these analytic functions. One 
such representation, with superior convergence properties, is 
a continued fraction ( C'-fraction ) which is approximated by 
truncation. The relation between these two representations 
is crucial for understanding of Pade approximants and is 
described below ED, 

For a given power series V(s) = Co + c\S + Cis 1 + ..., 
assume the existence of the reciprocal relation between the 
original series, as modified below, and a new series indexed 
by superscript V), 
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Now the original power series can be expressed as, 
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Next assume the existence of another reciprocal relation 
between the modified series from the denominator of the 
fraction in CD and a new series indexed by superscript ( 2 \ 
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This allows the expansion of the denominator of ifTTb in 
terms of another fraction, 
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By successively forming the reciprocal series we obtain a 
C-fraction, written in a compact form as. 
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By truncating the C-fraction in llTdl i we obtain its conver- 
gents which are rational fractions in s. For example the first 
4 convergents of (ITdl i are given as, 
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where c^ = — c 2 /c\ and c^ 2) = (c| — ciC3)/(ciC2). 


















The diagonal Pade approximant of degree M of V'(.s), 
hereafter appearing frequently in the text, is the (2M+l)th 
convergent of its C'-fraction representation in llTdl i. 
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In general a given analytic function can be approximated 
by PA [L/M\(s) where L and M are not necessarily equal, 
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Setting b 0 = 1, the denominator coefficients ?q , ..., b\,f are 
obtained by cross-multiplying (fT71 >. equating the coefficients 
of s L+1 ,s L+2 ,...,s L+AI to zero and solving the resulting lin¬ 
ear system. Next the numerator coefficients ao,ai,..aq are 
obtained similarly by equating the coefficients of s°,s 1 ,...,s i . 

Now consider the network of Figure Q] where the per-unit 
values of parameters in C3 are shown. The Taylor series 
for the unknown states, Vj (s) and V 2 (s) are obtained based 
on (O which are then used to compute the Pade coefficients. 
The concentration of zeros and poles of the diagonal Pade 
approximant, shown in Figure [2j defines the closest common 
branch point of Vj (s) and V^s) at Sb = 1.5 which is also 
given by Fabry’s theorem M as Sb — limn-j-oo Cn/Cn-\- 1* 
Here, as it is often the case for the class of problems in ([3]), 
analytic continuation by rational approximation is unnecessary 
as the power series already converge at s = 1 and thus are 
sufficient for obtaining Vj and V 2 . However, PA is a more 
efficient method as it converges to a given function at a much 
higher rate than the original power series does m. it can 
also discover the analytic structure of a given multi-valued 
function m. 

Bus 2 Bus 1 



III. Embedding the Voltage Magnitude 
Constraints in the Complex Plane 

For a generator node i G G C J\f, the real (active) power 
output. Pi = Re(Sj), is fixed whereas the imaginary (reactive) 
power, Qi = Im(Sj), is a free parameter which is adjusted so 
as to fix the magnitude of the voltage phasor Vi at a given 
setpoint value M$. Notice here the magnitude of a holomorphic 
function, Vj(s), is to be held fixed which forces it to be 
constant by open mapping theorem lU3l as the image of Vj(.s) 
in the complex plane is a subset of a circle and thus Vj can no 
longer be an open map. To resolve this contradiction we define 
an analytic function Vj(s) = ^^=0 s n independently of 
Vi{s) for i G G in such way as Vj(s)Vj(s) = Mf. Note that 
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Fig. 2: Zero-pole distribution of PA[100/100] 


7^ (Vi(s*))* and this distinction is the essential concept 
behind embedding voltage constraints and allows Vi(s) to 
adopt the value of V r at s = 0 with its magnitude approaching 
Mi as s increases. At s = 1, Vi(s) = (I^(s*))* for i G Q. 

Since for generator nodes (PV buses), Qi = Im(S'i) is 
unknown we define Si(s) and Si(s) as functions of complex 
variable s where S)(s) + Si(s) = 2 Pi and introduce 
as a complimentary set of equations to ©. At s = 1, © 
combined with (IT8l) sufficiently determine the AC power flow 
relations in a network with load and generation. Notice the 
different embedding of (1 1 Sab and (118bl) . The former ensures 
that Vj(s) has the trivial solution V r at s = 0. The latter 
enforces V t (s) = (Vi(s*))* at s = 1. Similarly Si(s) = 
3" s " 7 ^ (&i( s *))* for i G G except in their analytically 
continued form at s = 1. 

'«»>«* Weg(lSa) 
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Vi{s)Vi{s)=M 2 Vi GO (18c) 

Si{s)+ Si(s) = 2Pi ViGG (18d) 

The combined system of (|7]l and (IT8l) is adequately de¬ 
scribed by the following set of power series relations. 
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(a) Zero-pole concentration forms the Stahl’s compact set highlighting the common branch points of V\ (s), V± (s), Si (s) and V 2 (s) (Pi = 2.00). 
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(b) Transformation of the Stahl’s compact set on the feasibility boundary (Pi = 8.36) 

Fig. 3: Zero-pole distribution of PA[1000/1000] depicting the analytic structure of V\ (Y) (corresponding to Fig. 0. 


The key coefficients Cn (Vi € N— {r}) are progressively 
obtained by differentiating (1 1 9ab and (1 1 9bb with respect to s, 
evaluating at s = 0 and solving the linear system of (120ab and 
(120bb which itself requires the prior knowledge of dm. Cm, gm 
and g'm for m = 1 — 1. These coefficients are already 

obtained at previous stages through (120cl) - (120el ). 

Notice that = 2 Pi—g® whereas for rri > 1, gm = —gm. 
At s = 0, we obtain the trivial solution = V r for all nodes 
and subsequently c Y = Mf /cq ■ for generator nodes. 
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Fig. 4: 3-bus network with voltage control at bus 1 
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Now consider the modified network of Figure Q] where 
bus 1 has a generator that regulates its voltage magnitude at 
1.10 and generates Pi = 2.00. Figure [3a] shows the zero- 
pole distribution of the diagonal Pade approximant for V) (,s) 
forming the Stahl’s compact set E). The concentration of 
zeros and poles highlights the branch points of Vj (s) which, 
in this case, are common with Li(s), V 2 (s) and Si(s). From 
each branch point an analytic arc emanates and culminates 
in a different branch point or in a Chebotarev’s point of 
the Stahl’s compact set. The region of convergence is a disk 
bounded by the closest branch point Sb ~ —0.50. In contrast 
to the previous case, here, the concept of analytic contin¬ 
uation by Pade approximants is elegantly illustrated. Since 
linin^oo |c n+ i|/|c„| « 2, the coefficients tend to explode 
rapidly. Without Pade approximants based on these otherwise 
useless coefficients, it is impossible to recover the network 
solution. 

Figure [3b] shows the transformation of the Stahl’s compact 
set as Pi reaches the feasibility boundary. The branch point on 
the positive real axis has now moved to s = 1. Since past the 
branch point, there is a non-trivial monodromy, examining the 
PA solutions, as the degree of the diagonal Pade approximants 
is increased, reveals whether the power flow problem has 
a stable/high voltage solution or not. The location of this 
branch point can also serve as a proximity index to the 
feasibility boundary where the saddle-node bifurcation, i.e. 
loss of structural stability, occurs. 


IV. Superiority of PA over Newton-Raphson and 
Semidefinite Relaxation Methods 

Figure [5] shows a 7-bus network with 4 load (PQ) buses 
labeled 1-4, two generator (PV) buses labeled 5 and 6 and a 
reference (slack) bus. All values are in per unit. Line and load 
parameters are indicated as complex quantities. The generator 
voltage magnitudes are controlled at 1.10 and their active 
power output is 1.00. Newton-Raphson fails to solve this 
problem as it does not converge with a flat start, i.e. when 
initialized with all phase angles set to zero and all PQ voltage 
magnitudes set to 1.00. The first-order semidefinite relaxation 
also fails as it is not tight enough and the second-order 
(moment) relaxation finds an unstable/low voltage solution. 
In contrast PA method finds the desirable solution and, as 
the zero-pole distribution in Figure [6] clearly demonstrates, 
the operating point is on the stable branch and still has some 
margin to power flow infeasibility. 

Now consider the power flow in its polar form Q. The 
variables are phase angles of the 6 buses and voltage magni¬ 
tudes of the 4 load buses. The power flow has 10 equations, 
relating the active power of the 6 buses and the reactive power 
of the 4 load buses to phases angles and voltage magnitudes. 
To better contrast these methods we free a single parameter, 
Pq, the active power generated at bus 6. Each equation defines 


TABLE I: Pq =0.20 (all methods finding the stable solution) 


Voltage 

Pade 

Newton 

SDP 

SDP 

Magnitude 

Approx. 

Raphson 

(1st order) 

(2nd order) 

IVil 

0.9408 

0.9408 

0.9408 

0.9408 

IV2I 

0.9774 

0.9774 

0.9774 

0.9774 

\V 3 \ 

0.9953 

0.9953 

0.9953 

0.9953 

\v*\ 

0.9447 

0.9447 

0.9447 

0.9447 

TABLE II: 

Pq = 0.30 (failure of first-order relaxation 

Voltage 

Pade 

Newton 

SDP 

SDP 

Magnitude 

Approx. 

Raphson 

(1st order) 

(2nd order) 

|VT| 

0.9217 

0.9217 

- 

0.9217 

\Va\ 

0.9640 

0.9640 

- 

0.9640 

|Va| 

0.9897 

0.9897 

- 

0.9897 

|Vi| 

0.9403 

0.9403 

- 

0.9403 


TABLE III: Pq = 0.75 (false solution of moment relaxation) 


Voltage 

Pade 

Newton 

SDP 

SDP 

Magnitude 

Approx. 

Raphson 

(1st order) 

(2nd order) 

|Vi| 

0.7613 

0.7613 

- 

0.8504 

IV 2 I 

0.8658 

0.8658 

- 

0.9456 

|Va| 

0.9210 

0.9210 

- 

0.7960 

IV 4 I 

0.8888 

0.8888 

- 

0.1321 


TABLE IV: Pq = 1.00 (non-convergence of Newton-Raphson) 


Voltage 

Pade 

Newton 

SDP 

SDP 

Magnitude 

Approx. 

Raphson 

(1st order) 

(2nd order) 

|VI| 

0.5657 

- 

- 

0.8609 

IV2I 

0.7546 

- 

- 

0.9405 

|Vs| 

0.8394 

- 

- 

0.8178 

|Vi| 

0.8319 

- 

- 

0.1294 


TABLE V: Pq = 1.02 (false solution of Newton-Raphson) 


Voltage 

Magnitude 

Pade 

Approx. 

Newton 

Raphson 

SDP 

(1st order) 

SDP 

(2nd order) 

|Vi| 

0.5355 

0.1520 

- 

0.8575 

IV 2 I 

0.7380 

0.0673 

- 

0.9380 

\v 3 \ 

0.8283 

0.7376 

- 

0.8167 

IV 4 I 

0.8247 

0.7757 

- 

0.1296 


TABLE VI: Pq = 1.12 (non-existence of a physical solution) 


Voltage 

Pade 

Newton 

SDP 

SDP 

Magnitude 

Approx. 

Raphson 

(1st order) 

(2nd order) 

|VT| 

- 

0.1242 

- 

0.8363 

IV 2 I 

- 

0.0680 

- 

0.9234 

|V&| 

- 

0.7224 

- 

0.8086 

|Vi| 

- 

0.7609 

- 

0.1308 


a hypersurface in K n where n = 6 + 4+1. The intersection of 
these hypersurfaces, once projected onto the joint space of the 
freed parameter and a given variable, yields a series of curves 
in ]R 2 . Figure[7]shows these curves in the (P 6 , |X+j|) space. The 
stable/high voltage operating points of the network of Figure [5] 
can only be realized on the segment that is highlighted in red. 
This segment is consistently found by PA method for all values 
of Pq G [—0.114 1.057]. However Newton-Raphson and 

semidefinite relaxation methods concurrently find the stable 
branch only on a small subset of this interval (Table[U). Beyond 
Pq = 0.204 the first-order relaxation fails (Table [III). Beyond 
Pq = 0.739 the second-order (moment) relaxation finds the 
false branches (Table m. These branches are highlighted in 
green in Figure [TJi. Newton-Raphson convergence becomes 
erratic beyond Pq = 0.973 ('fable II Vk As Figure [7}r shows it 
either does not converge for certain values of Pq (Table IIVI) or 
it converges to low-voltage, physically unrealizable and thus 
false operating points (Table [V]). Beyond Pq = 1.057 there 
is no physically meaningful solution and the PA method re- 


































































TABLE VII: Erratic convergence of Newton-Raphson at P 6 = 1.0416...1.0424 « 1.042 


Pe 

1.0416 

1.0417 

1.0418 

1.0419 

1.0420 

1.0421 

1.0422 

1.0423 

1.0424 

\Vi\ 

0.3189 

0.4696 

0.2889 

0.1438 

- 

0.1438 

0.4887 

0.3210 

0.1436 

|W 2 | 

0.6256 

0.0520 

0.6478 

0.0675 

- 

0.0675 

0.7126 

0.6266 

0.0675 

|Va| 

0.7667 

0.8261 

0.6513 

0.7342 

- 

0.7342 

0.8123 

0.7671 

0.7341 

|Vi| 

0.7902 

0.8221 

0.1472 

0.7725 

- 

0.7725 

0.8147 

0.7904 

0.7724 


turns no solution whereas both Newton-Raphson and moment 
relaxation find false solutions (Table IVII) . It should also be 
noted that Newton-Raphson is not robust even in finding low- 
voltage solutions. This is shown in Table [Vi] and highlighted 
in Figure [8] where small perturbations at P e - 1.042 results in 
Newton-Raphson finding different branches or not converging 
at all. For industrial applications power flow parameters are 
typically expressed in 2, 3 and rarely 4 significant digits. 
Hence the set of values (1.0416,..., 1.0424) can be rounded 
to 4 significant digits and represented as 1.042 but applying 
Newton-Raphson to this set yields five topologically distinct 
solutions as well as non-convergence. Figure [8] shows the 
detail of the region highlighted by a dashed green box in 
Figures [TJi and [7]i. This region contains the closest bifurcation 


at Pe = 1.057 and presents a clear visual contrast between the 
performance of PA and those of Newton-Raphson and moment 
relaxation. Notice that Newton-Raphson finds operating points, 
mostly false, on all 6 solution branches as shown in Figure [8] 
whereas moment relaxation consistently finds the false branch 
that Newton-Raphson rarely discovers. 

V. Conclusion and Ongoing Work 
In this paper we have presented a method of solving the 
algebraic equations of AC power flow with voltage magni¬ 
tude constraints. We started by parameterizing the system of 
equations in the complex plane. This, under the assumption 
of the algebraic independence of the equations, renders the 
complex variables, i.e. voltage and power phasors, algebraic 











































































(a) Second-order semidefinite (moment) relaxation solutions (green) versus PA solutions (red) 



(b) Newton-Raphson solutions (blue) versus PA solutions (red) 

Fig. 7: Superiority of the PA method over Newton-Raphson and SDP methods shown in the context of the network of Fig. [5] 














Fig. 8: PA solutions (red) versus moment relaxation solutions 
(green) and Newton-Raphson solutions (blue) 

in s, the complex parameter, and hence expressible as power 
series in terms of s. These power series are in fact Taylor 
expansions around a trivial solution which defines the germ 
of the physically realizable and stable solution branch. The 
radius of convergence of these series is determined by the 
closest branch point of the algebraic curves that are given 
as the reduced Groebner basis of the parameterized system 
of equations. We demonstrated the pivotal role of Pade ap- 
proximants in analytically continuing the germ beyond the 
region of convergence of the series up to a point of non¬ 
trivial monodromy. This enabled the recovery of the physical 
solution to voltage and power phasors from the divergent 
power series. Equally important was the analysis of the zero- 
pole distribution of the Pade approximants that reveals the 
analytic structure of the complex variables. We illustrated how 
the concentration of zeros and poles form the Stahl’s compact 
set and that the branch point highlighted by the accumulation 
of zero-poles on the positive real axis serves as a proximity 
index to power flow infeasibility or voltage collapse. Thus 
PA method can definitively determine whether a physically 
realizable and stable solution exists or not. 

Next we demonstrated the superiority of the PA method 
over the conventional Newton-Raphson and the recently de¬ 
veloped semidefinite relaxation methods in the context of a 
7-bus network. Newton-Raphson, by far the most prevalent 
method in power industry, can be highly unreliable as the 
feasibility boundary is approached. It can also converge to 
solutions that are physically unrealizable and should thus be 
discarded as false. On the other hand, semidefinite relaxation 
methods are superior to Newton-Raphson in the sense that 
they do not exhibit erratic convergence behavior. We are aware 
of cases where, unlike the 7-bus network examined here. 


moment relaxation performs substantially better than Newton- 
Raphson in finding the stable branch. However, in general SDP 
methods can also fail, i.e. the relaxation is not sufficiently 
tight, or may consistently find the false branches. While the 
performance of these optimization methods can potentially be 
improved by better choices of objective function or adding 
extra constraints, we argue that the tendency to find false 
branches is a fundamental limitation of semidefinite relaxation 
methods which have their origins in real algebraic geometry 
and therefore cannot utilize the elegant concepts of analytic 
continuation and monodromy to solve problems of a complex 
analytic nature. 

We are currently investigating the generic structure of the 
Stahl’s compact set for the parameterized algebraic equations 
of power flow and the effectiveness of PA near the feasibility 
boundary. We will disseminate our findings in subsequent 
publications. This will also include a more expanded treatment 
of network controllers primarily the on-load tap-changing and 
phase-shifting transformers. 
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